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We examine analytically the existence and stability of phase-locked states in a weakly hetero- 
geneous neuronal network. We consider a model of A'^ neurons with all-to-all synaptic coupling 
^ l' where the heterogeneity is in the firing frequency or intrinsic drive of the neurons. We consider 

■ both inhibitory and excitatory coupling. We derive the conditions under which stable phase-locking 

is possible. In homogeneous networks, many difi'erent periodic phase-locked states are possible. 
Their stability depends on the dynamics of the neuron and the coupling. For weak heterogeneity, 
the phase-locked states are perturbed from the homogeneous states and can remain stable if their 
' homogeneous conterparts are stable. For enough heterogeneity, phase-locked solutions either lose 

stability or are destroyed completely. We analyze the possible states the network can take when 



^ ^ ■ phase-locking is broken. 

^ PACS numbers: 87.10. -fe, 87.22. Jb, 87.22.As, 05.45. +b 



I. INTRODUCTION 



Phase-locked and synchronized neuronal firing is found throughout the brain and this coherent activity may play an 
important role in behavior and cognition [Q-^. Given their importance, it is of interest to understand the mechanisms 
^ , responsible for their generation. Here, we consider a network of heterogeneous neurons that are coupled synaptically 
and determine analytically the conditions under which their firing patterns will be phase-locked. Our results could 
' have implications for synchronous phenomena in other areas of biology and physics that are based on pulse-coupled 
g ; oscillators |-|. 

[ — ' Much of the previous analytical work on phase-locking and synchronization in pulse-coupled neuronal oscillators have 
focused on homogeneous networks j^-^]. However, any biophysically relevant neuronal network will likely include 
■ heterogeneous effects. It is known that heterogeneity can greatly reduce synchrony in networks of phase-coupled 
H ■ oscillators [plfi^pTl] . Coexistence of synchronous and asynchronous states have been shown to exist in large networks of 



, heterogeneous pulse-coupled oscillators ||l8|,|T9| . Recently it has been shown numerically that heterogeneous assemblies 
^ , of neurons can fire coherently. However, the coherent activity can be greatly reduced even when the heterogeneity 
Q \ is very small [ pOpl| ]. Here, we examine analytically the existence and stability of phase-locked solutions of weakly 
^ ' heterogeneous neuronal networks. We consider an all-to-all coupled network of size N where the heterogeneity is in 

' the intrinsic firing frequency of the component neurons. We examine both inhibitory and excitatory connections. 

Phase-locking in homogeneous synaptically coupled neuronal networks has been studied analytically with different 



formalisms [|10|- |12| , |22| . Here, we generalize the formalism developed by Gerstner et al. to examine the existence 
and stability of phase- locked solutions. We find the conditions for stability of homogeneous phase-locked solutions and 
examine what happens as we allow the heterogeneity to increase. Our formalism can analyze all periodic phase-locked 
states but we mainly focus on the synchronized or in-phase state, the anti-synchronized state, and the splay-phase or 
anti-phase state. We find that stable phase-locking depends crucially on the time course of the synaptic coupling and 
on the response properties of the neuron to synaptic inputs. 

Neurons have recently been classified in accordance with their phase response properties lO] . Type I neurons have 
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the property that their phase response curves are always positive, i.e. the next firing time of the neuron is always 
advanced in response to a positive (depolarizing) input. In contrast, for Type II neurons the next firing time is either 
delayed or advanced depending on when the input is received. It has also been shown that Type I neurons admit 
arbitrarily low frequencies while Type II neurons do not ||2^ . Our results confirm previous observations that for Type 
I neurons, excitatory coupling is generally desynchronizing while inhibitory coupling is synchronizing pO|-p^. For 
Type II neurons, the opposite can be true although not necessarily. In previous work, the requirements on the rise 
and decay time of the synaptic coupling for stable phase-locking have been separately emphasized . We find 

that depending on the nature of the synaptic coupling conditions on both are required. A new result for homogeneous 
neurons is that the splay-phase or anti-phase state can be stable for inhibitory coupling for any finite sized network. 
For infinite TV, the splay-phase state is known to be unstable [p6|j29| ]. 

With weak heterogeneity, phase-locked solutions can exist but with small phase dispersion around their homogeneous 
counterparts. Stability is possible if the original homogeneous solutions are stable. The addition of heterogeneity 
cannot stabilize an unstable homogeneous phase-locked state. As the heterogeneity is increased, phase-locked states 
can either lose stability or cease to exist. Heterogeneity can also foster clustering. When stable phase-locking is lost, 
the network can enter asynchronous, harmonically locked or suppressed states. Fig. |l| shows examples of various states 
for two coupled inhibitory neurons obtained from numerical simulations of a biophysically based neuron model given 
in Appendix ^ 

In Sec. H we introduce the model and formalism we use for our analysis. In Sec. ttll we give the conditions which 
must be satisfied for a periodic phase-locked solution. In Sec. |^ we give general conditions for stability of locked 
solutions. In Sec. ^ we apply our conditions for the existence and stability of phase-locked states to two coupled 
homogeneous and heterogeneous neurons. We also examine suppression and harmonic locking. In Sec. VI we examine 



phase-locking and other dynamics in a network of N neurons. We discuss our results and conclude in Sec. VII 



II. MODEL 

Our network consists of N all-to-all coupled neurons. Each time a given neuron fires, a synaptic pulse is transmitted 
to all the other neurons. The effect of the synapse remains over a finite amount of time. The dynamics of neurons are 
generally described by biophysically based (Hodgkin-Huxley type) equations which are a set of differential equations 
for the current balance of the membrane potential and for the dynamics of the component ion channels [^^ . In 
Appendix ^, we give a set of biophysically based equations . We use these equations in numerical simulations to 
confirm our analytical results. These neurons are Type I in contrast to classic Hodgkin-Huxley equations which are 
Type II 

Generally, the differential equations for neuron models are complicated and cumbersome for theoretical analysis. 
Instead, we use the spike response method for our neuronal dynamics This formalism is fairly general and 

can be applied to biophysical neuron models. In this method, the neuron is considered to be a device which fires 
when a single scalar variable (membrane potential) crosses a threshold from below. The dynamics of the membrane 
potential consists of a sum of two sets of kernels akin to an integral expansion: 

I 3,1' 

where t\ marks the firing times of neuron i (i.e. the lih time neuron i reaches the threshold Oi) and j is summed from 
1 to excluding i. The kernels are nonzero only for t — tf > 0. The kernel rii{t) describes the intrinsic dynamics of 
the neuron after it has been triggered to fire. It represents the ensuing spike (action potential) and recovery behavior. 
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The kernel eij{t) represents the response due to synaptic inputs. We do not include the effects of self-coupling in (2.1) 
but can do so by modifying rii{t). Figure || gives examples of kernels. 
To linear order we can write eij (t) in the form 

e,j{t) = f Gis)hj{t~s)ds, (2.2) 
Jo 

where lijit) is the synaptic input to neuron i from neuron j, and G(t) is the linear response function. It has 
been shown that the linear response is often an adequate approximation for the full nonlinear response | |l^ , ^ , p7t . 
In biophysically based neurons, the synaptic input current is usually given by Iij{t) — gsSij{t){vi — Vs) where Vs 
is a constant reversal potential, gs is the maximal synaptic conductance, and Sij (t) is the synaptic gating variable 
obtained from a differential equation (see Eq. ( A4) in Appendix^. Models often use the simplification (t) ~ gSij{t) 



since (v.; — Vg) is approximately constant away from a spike The dynamics of the synaptic gating variable has 
previously been approximated by a recursive scheme where Sij (t) is updated each time the neuron receives an input 
with Sij(t) Sij{t) + Sf{t — P), where P is the time of the pre-synaptic spike, and Sf{t) is a stereotypical synaptic 
response [ p^pT|j29| ]. We call this type of synaptic model nonsaturating since Sij{t) increases without bound as the 



firing rate increases. On the other hand the synaptic model given by (A4) is saturating. Here Sij{t) saturates to a 



finite value as the firing rate of the pre-synaptic neuron increases. In the limit of very fast synaptic rise times we can 



use the approximation S{t) Sf{t — P) to model the dynamics of (A4) In this case the synapse has no memory 
of previous activity. 

In our network, we allow for heterogeneity in the thresholds of the neurons which implies heterogeneity in the 
intrinsic firing frequencies. We maintain homogeneity in the time courses of the kernels. We can make the threshold 
heterogeneity explicit with the transformation Vi —>■ Vi — li and 9i —>■ 9 — Ii, where li can be thought of as a normalized 
applied driving current, and is a constant threshold. We make the synaptic strengths of the neurons explicit by 
setting tij{t) — Je{t), where J represents the coupling strength. Heterogeneous coupling could be accounted for by 
setting J — Jij. When J > we call this excitatory coupling and for J < we call this inhibitory coupling. Putting 



this in Eq. (2.1) yields 



I 0,1' 

In this form, the network heterogeneity is explicitly expressed in the applied current li. We can set the threshold for 



all the oscillators to 1 by a simple rescaling without any loss of generality. In this formulation, Iij{t) in Eq. {2/1) can 
be thought of as the post-synaptic current while eij{t) can be thought of as the post-synaptic potential. Often we 
will refer to the rise time of the synaptic current which will mean the rise time of Iij{t). We note that the synaptic 
kernel eij(t) will itself have a rise time which can depend on both the rise and decay time of the synaptic current. 

Remark 1 We have assumed that the synaptic kernel does not depend on the time of arrival or phase of the synaptic 
input. For the integrate-and-fire model, as will be shown below, this is true. However, in general, there will be a 
dependence which we can express as e,;j = e{t ^t\^t — P„-^). For the ensuing analysis, we examine phase-locked solutions 
where the time since the last spike is fixed. Our assumption is then valid if |9„e(7i,w)| << |9„e(u,w)|. For neurons 
where the response changes drastically with phase, our analysis will need to be modified. 



1. Integrate-and-fire Model 

To compare to previous analytical results and to provide a simple concrete example we consider a leaky integrate- 
and-fire model |l0[-fi2|,p0[ with dynamics given by 
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^ - - V, ~Y.^{t- tl) + J2 JSfit - tj,), (2.4) 

I 3,1' 

where the firing times t\ indicate when Vi reaches the threshold value which we have set arbitrarily to 1. Each time the 
potential reaches threshold a delta function pulse is subtracted to reset the potential. As an example for a saturating 
synaptic current, consider a double exponential function 

^ f exp(-/3i) - exp(-ai), < i < t/, 

1 0, tf<t, 

where a > (3 and i/ is the time elapsed to the next synaptic input. For nonsaturating synapses, we can set tf — 
oo 0,0. 

Equation ( ^.4|) can be integrated directly flj,^ to obtain the kernels : 

= -exp(-t), t > 0, (2.6) 

, i4fl[exp(-/3i) - cxp(-t)] - j4^[exp(-ai) - cxp(-t)], 0<t<tf, 



[exp(-/3t/) - exp(-i/)] - [exp(-atf) - exp(-t/)] j exp(-(t - tf)), tf < t 

We show examples of the kernels in Fig. ^. Note that for a > P, e{t) > so that the time of the next firing is always 
advanced in response to a positive synaptic input and thus the integrate-and-fire model can be classified as Type I. 



III. EXISTENCE OF PHASE-LOCKED SOLUTIONS 



We first examine the existence of periodic phase-locked solutions. We derive a condition which all periodic phase- 
locked solutions must satisfy. Our condition generalizes that of Gerstner et al. [|l2| to include heterogeneity in the 
applied current and locking at arbitrary phases. Van Vreeswijk et al. derived a locking condition for the specific 
case of two homogeneous integrate-and-fire neurons locked at arbitrary phases. 

Our strategy is to assume that there is a periodic phase-locked solution and then check for self-consistency, i.e. the 
membrane potential must be 1 (the threshold value) at the time of the next firing. Consider each neuron labeled by 
i to be firing periodically in the past at t\ — {(pi — l)T, where I — 1,2, . . ., and < ipi < 1 is a. phase shift. Inserting 
these firing times into Eq. ( |2.3| ) gives 

v,{t)^I,+J2v(.t~{<j)^~l)T)+ J2 •Mt-i^i-l)^), t<(j,,T (3.1) 

1>1 j=ii,l>>0 

This assumes that each neuron will fire once within a period. We will discuss the situation where this does not occur 
later. 

Neuron i will fire next sX t — (piT. Inserting this into Eq. (^]^) and imposing self-consistency yields 

v,{(b,T) = l = I,+J2vilT)+ J2 J<l'T+{4>^-^j)T) (3.2) 

1>1 j=ii,l'>Q 

Note that for {(pi — (pj) < 0, the first term of the e kernel sum will be equivalently zero by definition of the kernels. 
The interpretation is that within the given period {I' = 0), neuron j spikes after neuron i and the synaptic effect does 
not appear until the next period, I' — 1. If Eq. (|]^) is satisfied then a phase-locked solution exists. The solution will 
be specified by the set of phases (pi and the period T. We note that the phase locking condition given by Eq. ( |3.2| ) 
consists of N equations for N phases and the period. However, the system is translationally invariant in time so one 
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of the phases can be arbitrarily fixed. This then leaves us with a well defined system of N equations and N unknowns 
( — 1 phases and the period T) . 

For a homogeneous system, the equations are symmetric with respect to the neuron index, implying a large number 
of possible solutions for the phases. The most notable solutions are the in-phase |]6|,p|-^ 16 or synchronized state 
where all the phases are the same and the splay-phase state pl[-p^ where the phases are spaced evenly apart over 
one period. There are also clustered solutions where equal numbers of neurons will fire at evenly spaced phases. An 
example is the anti-synchronous |lO|j3^ -|39[| state where half the oscillators fire a half a period apart. The inclusion 
of heterogeneity breaks the permutation symmetry. However, for small enough heterogeneity, near synchronous, 
anti-synchronous and splay-phase solutions that are perturbed from their homogeneous counterparts can exist. 



A. Phase-locked solutions for two neurons 



As a specific example, we consider the existence of periodic phase-locked solutions for the case of two coupled 
neurons. We consider neurons with heterogeneous applied current Ii > I2 and homogeneous coupling. We suppose 
the neurons have fired at times t\ = {(f)i — l)T, where ^ = 1, 2, 3, . . ., i = 1, 2. We choose the origin of time so that 
02 ^01, i-e. neuron 2 always fires after (or with) neuron 1. 



Applying condition (3.2) for A^ = 2 yields 

Viicl^iT) = l = h+ Y,WT) + MIT - ^T)], 
i>i 

V2{hT) = l=l2 + E + <^^)]' ^ 



(3.3) 
(3.4) 



/>0 



where 



61 > and we have taken into account that neuron 2 fires after neuron 1. We must solve Eqs. (3.3) 



and (3.4) for the phase (j) and period T. 



A relation for the phase can be derived by subtracting Eq. ( p.4| ) from (3.3) to give 

= /i - /2 - G(0) 



where 



(3.5) 



(3.6) 



i>i 



For homogeneous neurons, (/i = ^2), condition ( p.5[ ) for the phase is equivalent to that derived by Van Vreeswijk et 
al. ijl^. For e(0) = (which is generally true), we find that G'(O) — G(0.5) = G(l) = 0. In addition, G((/)) is symmetric 
around G(0.5). This implies that for homogeneous neurons = and (f) = 0.5 are phase-locked solutions which we refer 
to as the synchronized and anti-synchronized solutions respectively. If additional phase-locked solutions exist they 
must always appear symmetrically in pairs on either side of = 0.5. From Eq. (3.5), it is apparent that neurons with 
heterogeneous h cannot lock in a purely synchronous or anti-synchronous state. However, if the heterogeneity is not 
too large, there can still be phase-locked solutions near to synchrony and near to anti-synchrony. For heterogeneous 
neurons it is also possible that phase-locked solutions do not exist at all. As a specific example, consider the integrate- 
and-fire model. Using e{t) without saturation from (2.7) in condition (|3.6D, we find 



G(0) = 



J 



a 



„a<f>T-aT 



-a4>T 



-4,T' 



-aT 



J 



1-13 



J34>T-l3T 



a't'T-T 



-4,T 



-I3T 



(3.7) 



Figure H shows some examples of G(0) for various, T, a and [3. 
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To specify the period we add Eq. (|3.4[) to (|3.3[) and divide by 2 to obtain 



1 = / + Y^MIT) + -^HIT - 0r) + eilT -T + ^T)]}, 



i>i 



(3.8) 



where / = (/i + /2)/2. Using Eq. ( |3.§| ) and assuming the rise time of the synaptic current is very short compared to 
the decay time {a » [3), we find that the period for the synchronized state (0 — 0) for the integrate-and-fire model 
with saturating synapses satisfies [p8| 



f = /(I - e 



J{e 



-0T _ g-T\ 



(1-/5) 



(3.9) 



The period for nonsaturating synapses is given by Eq. (3.£) with an extra factor of (1 — e ^'^) in the denominator 
of the second term. The period does not change much for near synchrony where the phase is near zero. It has been 



found that (3.9) can be used to estimate the period of a synchronized inhibitory network of biophysically based Type 
I neuron models |28l. 



IV. STABILITY OF PHASE-LOCKING 

We examine the local stability of phase- locked solutions for a network of N all-to-all coupled neurons. Our formalism 
generalizes that of Rcf. [ p^ to include heterogeneity in the applied current li and for locking at arbitrary phases (pi. 
In this strategy, perturbations around the firing times of the phase-locked solutions are checked for growth or decay. 
Based on the perturbation dynamics we obtain sufficient conditions for stability. 

For neuron i we consider firing in the past at times 

tl^{<j,,-l)T + S,{n-l), 1 = 1,2,3,..., (4.1) 

where 5i{n — I) is a small perturbation around each firing time and n is the index for the current perturbation. We 
want to examine whether di{n) grows or decays as n — > oo. We will show that Si{n) depends on the perturbations 
of all of the previous firing times of all of the neurons in the network. To derive a mapping for Si (n) , we insert the 



perturbed firing times (}4.l|) into Eq. (2.3) to obtain 



v,it)=h + J2vit-i<j)^-l)T-^,in-l))+ J2 J<t~{4>o-l')T-Sj{n-l')). (4.2) 

The neuron will next fire at t = (piT + 5i{n). Hence we can set Vi{(j)iT + Si{n)) = 1, which implies after linearizing in 
di{n) that 

v,icP,T)^l-v,icj),T)S,{n), (4.3) 



where the dot denotes the first derivative with respect to time. We now set t = (f>iT in Eq. (4.2). Using Eq. ( |4.3| ) and 
linearizing in 5k{l) we obtain 

1 - v,{4>,T)5,{n) = I, +Y.^ri{lT) - 'n{lT)5,{n - /)] 

+ J[e{l'T+{^,~cl>j)T)~e{l'T+{cf>,-^,)T)6j{n~l')]. (4.4) 



For an unperturbed spiking history Eq. (|3.2[) holds. This allows us to simplify Eq. (4.4). Solving for 6i{n) yields 
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6i{n) 



Vi{4irT) 



E 

Z>1 



f){lT)5,{n ^ I) + 'J^i^'T^ 



cj,,)T)5,{n-l') 



where 



,T) = Y^fi{lT) 



i>i 



je[l'T + {(b^ - cj)j)T), 



(4.5) 



(4.6) 



Stability is determined from Eqs. (4.5) and ( [4.q ). It is not immediately transparent since a given perturbation of 
a neuron depends on all previous perturbations including those within the same cycle (i.e have the same index n). 
However, by turning the dynamics into a first return map we can obtain a sufficient condition for stability which we 
state in the following theorem. 



Theorem 1 // the series ) can be truncated then a sufficient condition for stability of a phase-locked state is that 



all of the coefficients of the series in (^) are positive 



Proof: We first make Si (n) depend only on previous firings by substituting for 6j (n) recursively on the right-hand side 
of the mapping (4.5) whenever (pi > (pj. This substitution can be done systematically if we index the neurons so that 
they are ordered by phase with (pi < (f)2 < ■ ■ ■ < Pn, i-c neuron 1 fires first, neuron 2 second, and so forth. With this 
phase ordering the mapping for 5i{n) is given by Eq. (4.5) with i — 1 and /' > 1. Since neuron 1 fires first it only 



depends on the previous firings of the other neurons. To obtain the mapping for S2{n) we must substitute for 6i{n) 
on the right hand side of mapping (4.5) with i = 2. We continue this recursive substitution scheme for all N neurons 



where for j = fc we must substitute for dj{n) on the right hand side of Eq. (4^) for j < k. We write the resulting 
mapping as 



Si{n) 



N 



l^^Ah,{l)5,{r 

j=i 1=1 



0, 



(4.7) 



where Mij{l) is obtained from applying the prescription described above and the sum over j now includes i. The 



elements of are composed of products and sums of the coefficients in the series (4.5). Thus, if the coefficients are 
positive then the elements of Mij are positive. Proposition 1 then directly follows from the following lemma. □ 

Lemma 1 // the series ^.S^ can he truncated then a sufficient condition for stability is given by Mij{l) > 0. 

Proof: By assumption the kernels ri{t) and e{t) vanish sufficiently quickly so that we can truncate the series in Eq. ( [4.5| ). 
This allows us to truncate the inner sum in Eq. ( [4.7D at some I = Im- We can then express the mapping ( [4.7|) in terms 
of a first return mapping by embedding in a higher dimension with 



where in block form. 



5{n) 



6{n) 
5{n- 1) 
5{n - 2) 



5{n) =M-5{n-l) 



( M(l) M(2) •■ 
I 
I 



(4.8) 



• m{iM-i) m(/m)\ 






(4.9) 



\5{n-lM + l) j V • • • I / 

I denotes the N x N identity matrix, d{k) is the column vector of length iV with elements 6i{k) and M(l) is the N x N 



matrix with elements Mij{l) from Eq. (4.7) 
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The dynamics of the perturbations are determined by repeated matrix multiphcations of JV[. The phase-locked 
solution is stable if all of the eigenvalues of Ai have absolute value less than unity. For simple enough systems we can 
calculate the eigenvalues explicitly (see Appendix |^) . This is usually not the case so we need to determine bounds 
on the modulus of the maximum eigenvalue p{M) (spectral radius). The spectral radius is less than or equal to any 
In particular p{M) < ||A^||oo = max^ ■ \Mij\. From Eqs. (|45|), (O) and (O) we find that 



given matrix norm | 

the rows of M(^) sum to unity, which implies that all of the rows M also sum to unity. Thus wc immediately find 
that 1 is an eigenvalue with eigenvector 5 — (1,1,1....,!). This corresponds to a global uniform shift in time and 
is due to the time translational invariance of the dynamics. It can be disregarded for stability considerations |l^. If 
the elements of A4 are non-negative then p{Ai) = ||A^||oo = 1- The spectral radius is the unique maximum modulus 
eigenvalue if A4 is a primitive matrix which is satisfied if A^™ > for some m |4^]. By inspection of A4 from ([l.S|), 
we see that M"^ > if Mij{l) > proving Lemma |^. □ 

It is important to note that Mij[l) > is not a necessary condition for stability. We make this explicit in the 
following corollary. 

Corollary 1 Stability is possible even if Mij has negative elements. 

Proof: If A4 is primitive we know that one eigenvalue is always 1, it is nondegenerate and all the other eigenvalues are 
within the unit disc. Stability is lost if an interior eigenvalue crosses the unit boundary, which can occur if we allow 
some of the elements of to become negative. Because the eigenvalues are continuous functions of the elements, 
we can make some of the elements negative and still keep the eigenvalues within the unit disc by the intermediate 
value theorem. □ 



Remark 2 Our perturbation scheme generalizes that of Gerstner et al. [g_2| who arrived at the stability equation ( |4.8| ) 
for the purely synchronized solution (all phases equal). They considered dynamics where e(0) = so only previous 
perturbations affect the current one. For a homogeneous network, they proved the necessary and sufficient condition 
for stability is Jj2i ^(^^) > 0- o^^' case with heterogeneity, pure synchrony is not possible. We are thus forced to 
consider locking at arbitrary phases and use the recursive substitution scheme above to calculate the mapping for the 
perturbations. 



A. Stability for two coupled neurons 



As an example we apply the formalism for N = 2. Equation (4.5) is 

Si{n) = ^ y[r,{lT)5i{n ~ I) + Je{lT - (j)T)S2in - I)], 

7)1 ^ ^ 



where 



Vl 



i>i 



S^in) = - y TlilT)62{n ~l) + — y e{l'T + ^T)5^{n - I'), 

f 2 V2 
^ 1>1 ^ /'>0 



(4.10) 
(4.11) 



^2 — 01. Since neuron 2 fires after neuron 1, we substitute Si(n) from Eq. ( I.IC ) into Eq. ( 4.11 ) to yield 

'i]{lT) .Pe{(t)T)e{lT - (t)T)' 



62in)^J2 



i>i 



V2 



ViV2 



52{n-l) 



Je[lT + (t>T) Je{(t)T)i){lT) 



V2 



VlV2 



6i{n-l) 



(4.12) 



Combining ( 4.1C ) with ( 4.12 ) we obtain 
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/ MIT) JeilT-<pT) \ 

■'^(0 ^ I ,Ji{lT+,j>T) Je{cl>T)ri(lT) i)(iT) J^ii)>T)eilT-<pT) I (4-13) 



which is to be used in Eqs. (^^) and (|4.9| ). The sufficient condition for stabihty of a phase- locked solution is given by 
Lemma ^ which we denote by M(^) > 0. 

For the case of two neurons we can also obtain the necessary condition for stability by examining whether or not 
one neuron can remain locked to the other one independently of whether they are firing periodically. This type of 
stability was considered by van Vreeswijk et al. [0. We derive their condition using the spike response formalism 
which we state in the following theorem. 

Theorem 2 A necessary condition for stability of a periodic phase-locked solution for two neurons is G'{4>) > 0, 
where G{(f)) is given by Eq. (S.t ) and prime denotes the first derivative with respect to (j). 

Proof: Suppose that neuron 1 fires at times t] — and neuron 2 fires at times tf — + 6* + 5{n — I), where 
I > 0. We consider the situation where neuron 2 has fired after neuron 1. Neuron f will satisfy the condition 
vi{ti) = 1 while neuron 2 satisfies the condition V2{ti + 9 + 6{n)) = 1, from which we derive the linearization 
V2{ti + 0) — 1 — V2{ti + 0)5{n). Inserting this into Eq. (2^) and linearizing in S{n — I) yields 

1 = /i + J2^viti - h-i) + Mti - ti^i ~9)^ Je{ti - ti^i - e)S{n - 1% (4.14) 
1 - V2{ti + 9)5{n) =I2 + - ti-i) - Viti - ti-i)S{n - I) + Je{ti - t2-i + 6)]. (4.15) 



Subtracting Eq. (4.15) from Eq. (4.14) and canceling terms, we obtain the condition for synchrony: 

/i - /2 = Jj2^eiti - t2-i + 6)- e(ii - ti-i - e)i (4.16) 



and the perturbation mapping: 



i>i 
where 



^^2(^1 +0) 



V2{ti +9)= Y^lvih - ti-i) + Je{ti - t2-i + 9)]. (4.18) 



Eq. (4.17) can be rewritten as 



where 



5{n) ^^M{l)5{n^l), (4.19) 



i>i 



M{1) = — [77(ti - ti-i) - Je{ti - ti-i - 9)]. (4.20) 



V2 

It will be shown in Lemma || that a necessary condition for stability is ^ M{1) < 1. 
Applying this condition leads to 



which implies 



< 1, (4.21) 
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J^[e(ii - h-i -9)+ e{h - t2-i + 0)] > 0, (4.22) 



For periodic firing we can set ti = {1 ~ l)T and 6 ~ (j)T in ( 4.22 ) which is equivalent to G'{(j>)/T > which leads to 
G'{(j)) > as stated in Theorem |. □ 

Lemma 2 For perturbed dynamics of the form 

(5(7i) = ^Af(0^(n-0, (4-23) 
1=1 

a necessary condition for stability is that '^M{1) < 1. If M{1) > 0, for all I, then the condition is sufficient as well. 



Proof: The dynamics can be expressed in matrix notation as in Eq. (4.8) where M(^) is replaced by the scalar M{1) and 
the identity matrices are replaced by 1. The corresponding matrix A4 is a companion matrix and the characteristic 
polynomial is given by p^ , pOt 

Im 

1 = ^M(0A-'. (4.24) 

1=1 

Let /(A) = J2i M{l)X~', so that eigenvalues of M satisfy /(A) = 1. Suppose J2i ^{1) > 1. Then /(A) > 1. However, 
/(A) ^ as A — > oo. Thus since / is a polynomial and hence continuous, /(A) = 1 at some A > 1, or in other words 
one eigenvalue will be greater than 1. Thus M{1) < 1 is necessary for all eigenvalues to be less than 1 and hence 
for stability. If M{1) > 0, then J^i -^(0 < 1 is necessary and sufficient for stability since /(A) < 1, for all A > 1. Thus 
there are no eigenvalues greater than or equal to 1. □ 

Remark 3 The condition G'{(j>) > is not a sufficient condition for periodic phase- locking. However, if the dynamics 
are constrained so that M{1) as given in Eq. ( 4.2C| ) is positive for all I, then J2i>i^{^) < 1 is the necessary and 



sufficient condition for arbitrary phase-locking. However, this constraint does not ensure periodicity. For instance, it 
may be possible that two neurons could be locked together in a higher order periodic or even an aperiodic orbit. 



V. TWO NEURON NETWORK 

We now apply our formalism to examine the existence and stability of specific phase- locked solutions for two coupled 
neurons. We will revisit the case of two homogeneous neurons and then consider neurons with heterogeneous applied 
current but with homogeneous inhibitory and excitatory coupling. In Appendix ^ we consider a simple example for 
two neurons where only the kernels from the last firing are kept. We compute the eigenvalues for stability explicitly 
and compare with our results for the general case. 

Applying the stability conditions for phase-locked solutions, we make the following proposition for the stability of 
two coupled neurons. Although we treat the spiking period T as a parameter in the proposition, it is important to 
note that it is an emergent time scale that depends on parameters of the network. 



Proposition 1 A periodic phase-locked solution for two coupled neurons is stable if m ( [^.J^ l i) i){lT) >Oforl>l, 



a) Je{4>T) > and Hi) Je{lT ± (f)T) > 0, for all I > 1. These are sufficient but not necessary conditions. 



Proposition]^ is established by analyzing the matrix elements of M(Z) in ( 4.13| ). By assumption vi > and V2 > 



Conditions i), ii) and iii) then imply that M(Z) > proving the proposition by Lemma Condition i) is generally 
true for most neuron models. Condition ii) is strongly dependent on the neuron type. Condition iii) is contingent 
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on the neuron type as well as the neuron parameters, i.e. strength of the applied current, the synaptic strength and 
the synaptic decay time. Gerstner et al. proved Proposition |l| for the homogeneous synchronous case (0 = 0) for 
N coupled neurons where N is large. We discuss this case in Sec. |VI A . Condition iii) corresponds to a special case 
of their locking theorem. This puts a condition on how slow the onset of the synapse can be. As will be discussed 
below, condition ii) was implicitly assumed in their theorem. Van Vreeswijk et al. [|lO| emphasized condition ii) in 
their criterion for stability. 



A. Homogeneous Neurons 



The existence and stability of phase-locking for two homogeneous neurons has been analyzed in detail by van 
Vreeswijk et al. and by Hansel et al. [0 for excitatory connections. These analytical results used phase reductions 
or considered integrate-and-fire dynamics. Our formalism is applicable to a more general class of neuron models. In 
Sec. VE , we consider heterogeneous neurons. Gerstner et al. considered the stability of the synchronous solution 
for a large homogeneous network. We will relate the stability conditions of Ref. [0| to those of Ref. |T3]. Periodic 



phase-locked solutions are given by Eqs. (3.5) and ( |3.8[ ) for Ii = I2. As discussed in Sec. IIIA, these include the 
synchronous {(j) — 0) ^^id anti-synchronous (0 = 0.5) solutions as well as a possible third locked solution. 



1. Inhibitory Coupling 

Here we relate the results and formalism of van Vreeswijk et al. [|l0| to our formalism. First consider synchronous 
locking. As mentioned, condition i) in Proposition |^ is satisfied for most neuron models. For Type I models, e(t) is 
always positive so condition ii) can only be satisfied if e(0) = 0. This condition is equivalent to a nonzero rise time of 
the synaptic current which was emphasized by Van Vreeswijk et al. [ p^ as the condition required for synchrony with 
inhibition. In order for synchrony to occur with inhibition, the effect of the first neuron on the second one cannot 
be felt immediately. They gave an empirical criterion that the synaptic rise time be slower than the width of the 
spike. For certain Type II models where e(i) is initially negative and then becomes positive Je(0) could be greater 
than zero (see Fig. ^). Condition iii) is easily satisfied for Type I neurons (which includes the integrate-and-fire 
model). Recall that inhibitory coupling is given by J < 0. Thus, if e[t) is positive with a maximum at t ~ tc and is 
monotone decreasing for t > tc (which applies to Type I neurons) then condition iii) in Proposition |l| is equivalent 
to T — <j>T > tc. The decay rate of e{t) is controlled to some extent by the decay rate of the synaptic current as well 



as the intrinsic membrane decay rate (see Eq. (2.7)). The period on the other hand is an emergent property of the 
network which depends on the applied current, synaptic strength, and synaptic time scales j2^. For Type I neurons 
it can be estimated from Eq. ( |3.9D . For certain Type II neurons it may be possible that condition iii) can never be 
satisfied so synchronous solutions would not be possible with inhibition. 

Proposition gives sufficient conditions for stability. By Corollary^ synchrony may be possible even if Proposition^ 
is violated. To find when locking is not possible we examine the necessary condition given by Theorem From 



Eq. (3.6) we note that the necessary condition for stability 

G'{0) = JTe{0) + 2JT^e{lT) > 0. (5.1) 



i>i 



is satisfied for Proposition]^ as expected. If e(0) > —2 J2i>i ^('^) then synchrony is unstable. Satisfying this condition 
depends on e(0) which is governed by the rise time of the synaptic current. In order for stability, the neuron must 
be able to fire without being immediately and strongly inhibited by the synapse. This relates to the van Vreeswijk et 
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al. statement that the rise of the synapse must be slower than the generation of a spike for stabiHty to be possible. 
For the integrate-and-fire model with nonsaturating synapses, a zero rise time guarantees instability. This can be seen 



by examining Eq. (3/7) where it can be shown that G"(0) < if a — >• oo and T > 0. This verifies previous results of 
local instability of the synchronous state for zero rise time with inhibition in the integrate-and-fire model . 
By Proposition |l|, the anti-synchronous solution is stable as long as the synapse is fast enough so that e{0.5T+lT) < 



0, for ^ > 0. From ( p.q ) we find the necessary condition for stable anti-synchrony is given by 

G'{0.b) = 2JTJ2KIT -0.5T) > 0. (5.2) 



i>i 



Condition (5_^) shows that G"(0.5T) could be negative if the onset of inhibition is very slow. For the integrate-and-fire 
model this requires that a and P be small enough in Eq. ( p.7[ ) (See Fig. ||). This is in accordance with the results of 
van Vreeswijk et al. | pO[ . 

The synchronized and anti-synchronized solutions can both be stable simultaneously if the rise time is nonzero and 
e{t) is monotone decreasing for t > 0.5T. If this is the case then additional phase-locked solutions must also exist. 
Since G{(j)) is continuous, if G'(0) is positive at 0, 0.5 and 1 then it must cross zero at least once between and 
0.5 and between 0.5 and 1. If it crosses only once then this solution must be unstable. (See Fig. ||). Since G{(j)) is 
symmetric about (j) = 0.5, we can think of these two solutions as one. This third solution was obtained numerically in 
the integrate-and-fire model with nonsaturating synapses in Refs. |l0|jl| where it was found to arise from a pitchfork 



bifurcation of the anti-synchronous solution. It approaches the synchronous solution as the synapse becomes faster 
compared to the period (a and /3 get larger) . The location of the third solution with respect to the synchronous and 
anti-synchronous demarcates the basin of attractions for the two stable solutions. All phase- locked solutions can be 
unstable if the onset of the synapse is so slow that e{t) is still increasing at t = T because M is no longer a positive 



matrix. For slow enough synapses, the necessary conditions for stability given by (5.1) and (|5.2|) can both be violated 



2. Excitatory Coupling 

The case of slow excitatory coupling for two neurons has been analyzed previously by van Vreeswijk et al. pC| ] 
and Hansel et al. pd| ]. They both show that for Type I neurons, synchrony is generally unstable and anti-synchrony 
is stable for slow synapses. Our analysis agrees with their results. For two neurons, stability is not ensured since 
e{lT) < 0, violating Proposition]^. From Eq. (5.1) we find that the necessary condition for synchrony in Theorem|| is 



violated for excitatory coupling unless e(0) is positive and large. An examination of G"(0.5) for J > in (5.2) shows 
that anti-synchrony is stable if the rise time of the synaptic current is very slow. 

However, stability might be possible for Type II neurons with excitation as pointed out by Hansel et al. [pl| . This 
can be seen in our formalism because e{lT) could be positive for Z > 1 if the phase response curve had both positive 
and negative parts. We show an example in Fig. |[ For the Hodgkin-Huxley equation, Kistler et al. have shown 
that the response of the membrane to a synaptic input can behave in a way to support synchrony. 

For the integrate-and-fire model with instantaneous (delta function) synaptic coupling, stable synchrony is possible. 
To see this consider the integrate-and-fire model (p]4) with a decay constant 7 for the membrane potential restored. 



This implies that ri[t) = — exp(— 7i) and e(t) = exp(— 7<). Thus e{t) = —jexp(—"ft), which implies Afi2 in (4.13) is 
negative. However, these negative elements can be made arbitrarily small if 7 and J are made small and stability is 
possible by Corollary |[ For a combination of small coupling and weak decay, the synchronous state can be stable. 
This is in accordance with the result of Peskin ||3^ who showed that the product J7 must be small for synchrony in 
the integrate-and-fire model. 
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For the integrate-and-fire model, the third phase-locked solution which comes from a pitchfork bifurcation of the 
anti-synchronous solution is always stable ■ As the synapse becomes faster and more pulse-like the third solution 

approaches the synchronous solution. Thus for very brief synapses, two neurons locked in the third solution, can be 
considered to be in a state of almost synchrony Q . 

Remark 4 It is important to note that our local stability is different from the global stability of MiroUo and Stro- 
gatz 11] and related work for pulse-coupled oscillators. In these papers, an absorption condition is assumed for 

synchrony. By this it is meant that if two oscillators fire together within a certain time window, then they are consid- 
ered synchronized. These papers show that for pulsatile excitation and almost all initial conditions, the oscillators will 
eventually get close enough to be absorbed. What we have examined here, which is in line with Refs. [|0 
local stability of coupled oscillators that still influence each other even if completely synchronized. After the oscillators 
have been absorbed they may not synchronize in the local sense. They could be in a state of almost synchrony, make 
small oscillations around synchrony, or even be repelled away |l2|j4l| ]. 



B. Heterogeneous neurons 

We now apply our formalism to two coupled heterogeneous neurons. In this case, phase- locked solutions need not 



exist. For large enough heterogeneity or weak enough synaptic strength, condition (3.5) may not be satisfied (See 
Fig. ^). When there are no phase-locked solutions, several possibilities can arise. If no stable phase-locking exists and 
both neurons are able to fire within a period then a possible outcome is a state of asynchrony, where each neuron is 
effectively decoupled and fires independently of the other neuron. For strong inhibitory coupling, we can get a state 
of harmonic locking where the neurons are locked at some ratio of their periods. The locking period can be very long. 
If the inhibition is so strong so that the slower neuron never fires then we call this state suppression |2|. SeeFigj^for 
numerical examples of these states. We will discuss harmonic locking and suppression in more detail in the following 
sections. 

Under certain conditions, stable phase-locked solutions can exist. For weakly heterogeneous neurons these solutions 
will be near to the homogeneous locked solutions. We concentrate on the near synchronized and near anti-synchronized 
solutions. From ( |3.5[ ) we see that for a given level of heterogeneity, | J| must be large enough for a phase- locked 
solution to exist (See Fig. H). The phase (p is determined by the heterogeneity and the coupling strength. For small 
heterogeneity and small phase, (f) varies linearly with Ii — I2 and inversely with ~J. 

For excitatory coupling, the stability results for homogeneous neurons can be applied to the heterogeneous case. 



Near synchrony is generally unstable for Type I neurons. For integrate-and-fire neurons, by analyzing M(Z) in (4.13), 
we see that stability may be possible in the case of weak coupling, weak decay of ri(t), and extremely brief synapses 
as in the homogeneous case. The conditions for stability of the near anti-synchronous solution is possible if the 
heterogeneity is mild and the synapse is slow. Near synchrony can be stable for heterogeneous Type II neurons by 
Proposition 

For inhibitory coupling ( J < 0), stability of near anti-synchrony can be obtained directly from Proposition |l| but with 
(j) shifted slightly from 0.5. As long as conditions ii) and iii) arc still satisfied then anti-synchrony is stable. Stability 
of the near synchronous case where is small but nonzero will depend on the neuron properties. If Je(0r) > 0, then 
stability is ensured by Proposition |l|. For Type I neurons this will not be generally true but e(0T) — may hold as 
a result of delays in the onset of the synapse or from refractoriness of the neuron. (We note that Je{(j)T) < for 



the integrate-and-fire model (2.4).) For Je{4)T) < 0, stability is not guaranteed because elements of the second row 



of M(Z) could become negative. However, by Corollary |l|, negative elements do not immediately imply instability 
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because the eigenvalues can remain within the unit circle even if M(/) takes on negative values. This leads to the 
foUowing proposition on the stabihty of inhibitory heterogeneous neurons. 

Proposition 2 Stability of near synchrony for Type I inhibitory coupling is possible if the heterogeneity and strength 
of the applied currents are both weak enough, and the synaptic strength is strong enough. 



We show this by examining the elements of M(^) in ( 4.13 ) individually and then determine conditions for them to 



be positive. As before, we assume that rjilT) > 0, for all I. We also assume that e{t) has a single maximum at t — tc 
and that T - (t)T > t^ so that e{lT ± (j)T) < 0, for all I. This impHes that Mii{l) > and A/i2(0 > 0. Heterogeneous 
phase-locking implies that (f> > which affects elements M2i{l) and M22(0- 
First consider 

M2,(0^^^^^^^^ + «^. (5.3) 

For J < 0, this element is positive if the magnitude of e{4>T) > is small enough and that of e(/T + 0T) < is large 
enough and it decays away as slow or slower than fj{lT). If we assume that e(0) = then £{<j>T) small implies that 
(j)T must be small. But if T is too small then e{lT + 0T) can become too small or even positive. The phase can be 
made smaller if the level of heterogeneity is reduced or if the synaptic strength \J\ is increased. From Eq. ( |3.9| ) we 
see that the period T is made longer if the synaptic strength is increased or if the applied current is reduced. Thus 
in order for Af2i(0 to be positive we require weak enough heterogeneity and strength of the applied currents and a 
strong enough synapse. 
Next consider 

^^^^^^^^^JMWOT^, (5.4) 

V2 V1V2 

This element is positive if e{(f)T) > is small and the decay of e{lT — (j)T) < is as fast or faster than fj{lT). Thus 
requiring M2i(^) > and M22{1) > implies that e{lT ± (j)T) and ri{lT) must have the same decay rates for / > 1 
or there will be negative elements in the second row of M(Z) for some I. For saturating synapses and fast rise times, 



this can be satisfied (see Eq. (2.7)). The decay rates will in general be different for nonsaturating synapses. However, 
even if this is the case, as long as the decay rates are fairly close or fairly fast, the negative elements can be made 
arbitrarily small and stability is still possible by Corollary |l|. This result indicates that synchronization should be 
easier to achieve with saturating synapses as opposed to non-saturating synapses for heterogeneous neurons. 

Remark 5 Proposition || confirms the numerical observations in Ref. where it was found that synchrony was 
only observed in what we called the phasic regime. It was shown in Ref. |Q that the conditions of Proposition || lead 
to the phasic regime where PT > 1. We can use this to estimate the boundary between synchrony and asynchrony 
in the J - (3 parameter plane. From our above analysis and previous numerical simulations pl[ | we find that /?T is 
approximately constant in the synchronous regime. At the boundary, for the level of heterogeneities chosen (5% in 
intrinsic frequencies), we found that /3T ^ 1 defined the boundary. Recall that the period of a synchronized network 



with saturating synapses approximately satisfies the transcendental relation (3.9) [p8|, which we rewrite as 



(e ^ — e ^) 



where we have substituted T ~ r = /3 ^ . Eq. (5^) gives the synchrony- asynchrony boundary for g as a function of 



(3. We can obtain a simpler expression if we let that t 1 + 5 and expand for small 5: 

J - -[/(e - 1) - e + /(5.] (5.6) 
Thus we see that J depends linearly on 5 r — 1. This was observed in numerical simulations [ pT[ . 
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C. Suppressed States 



In the previous analysis for phase-focking, it was assumed that there was a phase-locked solution with period T and 
phase (/). However, for inhibitory coupling this may not always be true. It may be that some neurons never fire at 
all. We called such a state suppression in Ref. We illustrate this for two neurons. Suppose Ii > I2 and neuron 
1 fires before neuron 2. If neuron 2 never fires (i.e. limt^oo ^2(i) < 1) then we say neuron 2 is suppressed. As we will 
discuss in Sec. VI C , suppression can also occur in networks of N neurons. 

Here we derive the conditions for suppression for two inhibitory neurons. We assume that neuron 1 has been firing 
periodically in the past dXt={n — l)T, where / > 1 and nT denotes the current time. We also assume that neuron 2 
never fires. We let J = —g, where g > 0- Using Eqs. (2^) and ( |3.2| ), this can be expressed as 

ViinT) = 1 = h +y^r,ilT), (5.7) 

(5.8) 



lim V2{nT) = l2-gye{lT)<l. 



1=0 



Equations ( p3.7[) and (5^) give the conditions for suppression. 



Using the kernels ( p.6[) and (2/7) for integrate-and-fire neurons with saturating synapses and using the condition 

and obtain the condition 



a >> P we can compute the sum explicitly in Eq. (5 

9 



-(3T 



h 



1-/3 (1-e-^) 

We see that suppression can be achieved with large enough g 
Eq. (5.7) the period is given by 



< 1. 



(5.9) 



T = In 



h 



or small enough T. Without self-inhibition, from 



(5.10) 



If we include the effects of self-inhibition (i.e. the neuron inhibits itself with strength g when it fires), the period T is 
given by Eq. (3.9) with I = Ii. We can then use Eq. ( |3.9| ) to simplify Eq. (5.9) and obtain 

1 



1 



< 1 



(5.11) 



From this expression we see that if the neurons are homogeneous (/i = I2) then suppression is not possible. For 
heterogeneous cells, if I2 < Ii then suppression will definitely occur for large enough T. The longer the period the 
more likely suppression will occur as opposed to the case without self-inhibition. 

We can estimate the transition boundary to suppression 'm g ~ (3 space. Without self-inhibition, the period of 
neuron 1 is not affected by the synapse and we can assume that the boundary to suppression is given by condition 
( ^.9D with an equality 

-/3T _ „-T 
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1-/3 (1 



1 



1, 



(5.12) 



where the period is given by (5.10). 

With self- inhibition, the period now also depends on the synaptic parameters. From ( 5.11 ) we see that at a critical 
period suppression will occur. The boundary will then be given by the conditions specifying that period. Consider 
the case where where T >> 1 but j3T > 1. This corresponds to what we have called the phasic regime |^l|]. This 
allows us to ignore factors of e~'^ in the period condition (3.9). The suppression boundary is then given by 
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for fixed T such that condition (5.11) is minimally satisfied. Solving for g we then obtain 

g^(/i-l)(l-/3)e''^. (5.14) 
We see that for /? < 1 and T large enough, g increases as /3 increases as was observed numerically 



D. Harmonic locking 



It is also possible that two neurons may be harmonically locked with a period ratio pi : P2- The ratio need not 
be rational. Suppression can be considered to be harmonic locking at the ratio 1 : oo. Consider two neurons locked 
in the ratio of pi : p2 so that neuron i will fire pi times before repeating its firing pattern. Let the length of this 
periodic cycle be T, which is the same for both neurons. See Fig. |l| for an example of a 4 : 3 harmonically locked 
state. Consider two inhibitory neurons firing in the past at 

I 



— 

where [•] indicates least integer, I — 0,1,2,..., and 
firings. In the next cycle, neuron i will fire at t — T - 



T 



mod Pi ' 



(5.15) 



/ mod Pi 



is the phase advance per firing which repeats after pi 



where m,- takes the values 0, 1, 



1. If we substitute 



this into Eq. (2.3) we obtain 

V^{T- ■ " 



l>7ni 







(t + 









T 



I mod Pi 



l'>n 



T 



^f'modp, -'/'k)^ 

(5.16) 



where n is the index for which 



n mod pj 



> is minimal (i.e. n gives the index of the nearest firing of the other 
Pi + P2 equations that must be solved simultaneously for the period T 



neuron). Equation (5.16) is a system of M 
and the phases Because of the freedom to shift globally in time we can fix one of the phases (e.g. set = 0), 
which then leaves M equations for the remaining M — 1 phases and the period T. For a fixed set of parameters, there 
is an infinite set of possible harmonic states with different period ratios. 

Stability of the harmonic state can be established by applying results of the previous section. A stability analysis 
can be performed as before by perturbing the time of each firing and constructing a mapping for the perturbations. 
The mapping will map all the perturbations from one cycle at < = kT to the next cycle at i = (fc + 1)T, where k is 
the least integer oi l/pi. The sufficient condition for stability will be given by Theorem ^ We can examine the local 



stability of the relative firing times of the two neurons. As shown in Sec. IV A, for inhibitory coupling, a sufficient 
condition for stable phase-locking is that the phases between two consecutive spikes be either very small or far apart. 
More specifically terms such as gi{kT -\- [(f)\ ^^^^ ^, — (j)\n.)T) will arise in the perturbation series. The terms with fc = 



correspond to the term ge{(j)T) in Sec. VB . Recall that stability required gi{4>T) to be small. Terms for fc > will 
be positive if the synapse decays fast enough. This then implies that stable harmonic locking will involve spiking 
patterns where the spikes from the two neurons are generally either coincident, close together or fairly far apart (See 
Fig. ||. This was observed numerically |^ . 

Numerically, harmonic locking was more likely to be observed if self-inhibition was included pll . One reason may 



be that the locking conditions (5.16) are easier to satisfy with self- inhibition because the equations become more 
symmetric. Another reason may be that the transition from synchrony to suppression is slower as the parameters 
are changed for neuronal networks with self-inhibition compared to those without it, leading to a larger region in 
parameter space where harmonic locking is observed. 



16 



VI. iV COUPLED NEURONS 



We now examine the dynamics of a network of TV neurons. Many of the results from two coupled neurons can be 
carried over to the iV-neuron network. For N homogeneous neurons there are many possible phase-locked solutions. 
As discussed in Sec. [II, all symmetric combinations of the phases are solutions. Theorem |l| gives the condition for 
stability. We examine the existence and stability of synchronous, splay-phase, clustered, and unlocked states. Large 
homogeneous networks have been examined previously using different formalisms [ pl]p^jl5| ] . We consider the effects 
of weak heterogeneity. It is known that strong heterogeneity can induce various complicated states |lq,p^. 



A. Synchrony 

For a homogeneous network, the synchronous state is a phase-locked solution. Stability has been examined in the 
limit of large N by Gerstner et al. (See Remark ||). They proved that as long as the influence of the synaptic kernel 
is rising at the time of the next spike then synchrony is stable, i.e. X]i>i JK^T) > 0, where J < ( J > 0) corresponds 
to inhibitory (excitatory) coupling. Implicit in their proof is that the synaptic kernel does not act immediately (i.e. 
e(0) = 0). 

With the inclusion of heterogeneity, the phases are shifted apart for a locked solution. However, for weak hetero- 
geneity, a stable near synchronous solution is still possible. This was observed in numerical simulations ||2^. Here 
we illustrate that for inhibitory coupling, if the neurons are indexed so that Ii > I2 > ■ ■ ■ > In then there could be 
a stable phase- locked solution where the phases are arranged as 4)i < (t>2 < ■ ■ ■ < (t)N provided the heterogeneity is 
small and close to uniformly distributed. We translate time so that the phases are all positive and small. 

First consider a homogeneous system with a period satisfying the condition 

1 = ^" + Y.iim - 9{N - l)t{lT)]. (6.1) 

where J — —g, g > 0, and / = (1/A^) J2j ^^^^ ^'^^ Ii = I + ^i- Since the applied currents are all positive then 

Al > and Ajv < 0. If a phase-locked solution exists in the heterogeneous system it must satisfy the condition 

l^I + A, + J2vm- 9^{lT+{ct^^-4>,)T)-Y,ge{{<j,,-4>,)T), (6.2) 



Subtracting condition (6.1) from (|6.2|) gives 



A, = ^ g[e{lT+{<i>,-<l>j)T)-e{lT)]+Y,ge{{cj>,-<i>,)T) (6.3) 



Consider condition (6.3) for neuron 1 (i 



Ai= ^ g[e{lT+{ct>i-ct^j)T)-e{lT)]. (6.4) 



Since we have assumed that 0i < <f>j^ then condition (6.4) can be satisfied if e{t) is monotone decreasing for t > 



T + {(pi — 4)n)T- Now consider condition yjpM for neuron N 



An^ J2 gWT+{4>N~4>j)T)-e{lT)]+Y,ge{{4>N-4>o)T) (6.5) 

j<N,l>l j<N 



In this case Ajv < and (pN > <Pj- Condition (5.5) can be satisfied if e{t) is monotone decreasing for t > T 



3W 



(j)N-i)T and e{{4>N ~ 4'j)T) is not too large. We can perform a similar analysis for all the other neurons. The result is 
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that as long as e{t) is monotone decreasing for i > T + — 4>n)T and it does not rise up too quickly for small t then 
near synchrony is possible with the phases arranged as 4>i < (j)2 < ■ ■ ■ < 4>n, i-e. they will all fire in rapid succession. 
We should note that the period is controlled by the applied currents, the synaptic strengths, and the synaptic time 
course. 



Stability is determined from the perturbation series for stability (4.5) which will be composed of terms with coeffi- 
cients of the form r){lT), —ge{lT+ {(f>i — (f>j)T), and —ge{{(j)i — (t>j)T), for all i, j and I. By Theorem |l|, the phase-locked 
solution is stable if all of the coefficients are positive. The monotone decreasing condition guarantees that the first 
two sets of coefficients are positive. For a homogeneous system, the third set is zero. For a heterogeneous system, 
the third set can be negative for Type I neurons. However, by constructing the return mapping (4.8) it can be 
shown that stability is possible as long as they are small by using Corollary |l| and generalization Proposition ^ to the 
heterogeneous N neuron network. 

In the case of Type I excitatory coupling and zero rise time of the synaptic current, synchrony can be globally stable 
but it can be locally unstable If a neuron fires late with respect to the others then it will be pushed forward 

and eventually catch up to the synchronous population, but if it fires early it will be pushed away until it catches up 
to the next cycle. MiroUo and Strogatz |^ considered an absorption condition so that once two neurons fire within 
a given time window they are considered synchronized and remain synchronized. This absorption condition ensures 
global stability. 



B. Splay-phase and Clustered States 

For a homogeneous network, the splay-phase state also known as anti-phase, 'rotating wave' or the more colorful 
'ponies on a merry-go-round' [0,E5l|3^-Q is a possible phase-locked solution. In this state, all N neurons in the 



network fire with a phase separation of 1/N. This solution can be verified by examining condition (B.2). With weak 
heterogeneity, a near splay-phase state should exist by a similar argument to the perturbed synchronous case. 

We can show stability of the splay-phase state with an N neuron version of Proposition |l|. We note that the 



coefficients of (4.5) are composed of factors of the form r]{lT) and Je{lT ± jcjyT) where = 1/-/V, and j = I, . . . , N 
is the neuron index. If these factors are positive then there is stability by Theorem |l]. In the excitatory case, the 
splay-phase solution can be stable if the rise time of the synapse is very slow so that e{jcj)T) > 0. In this case the 
splay-phase is stable for similar reasons as to why the anti-synchronous solution is stable for two neurons. This has 
been noted previously |l^Jl^. With weak heterogeneity, stability of the near splay-phase state can still be stable 
provided the synapse is slow enough. 

The splay-phase state can also be stable for Type I inhibitory connections if the synapse has a fast enough rise time 
so that i{T/N) < and e{t) is monotone decreasing for t > T/N. Stability may be affected if the action potentials 
(spikes) of the individual neurons overlap in the splay-phase state. This can occur if the network is so large that 
T/N is on the order of the spike width. We have verified the existence of the splay-phase state numerically in our 
biophysical model for = 4. The stability of the splay-phase is interesting because in the mean field limit {N oo) 
it has been shown to be unstable for inhibitory coupling 26 4^. In the mean field limit for homogeneous neurons, 
the splay-phase state corresponds to an 'asynchronous' state where the density of the firing times of the neurons is 
constant. We can see how this instability arises in our formalism. As more neurons are added, the synapse must rise 
faster for stability. In the limit of infinite numbers of neurons, stability is not possible no matter how fast the synapse 
since {T/N) — > 0. With weak heterogeneity, a stable perturbed splay-phase state is possible as long as the rise time 
of the synapse is shorter than the smallest phase difference. Stability should be more affected by heterogeneity as the 
network size increases since the phase differences become smaller. 
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There is also the possibility of phase-locked clustered states. Clustering in neuronal networks has been investigated 
before (l|Jl|,||-||. For homo geneous networks, groups of neurons can fire together in unison but at different phases 
with other groups. Clustered states can be shown to be stable if e(0) = and Je{lT + {(j)a — (j)b)T) > 0, where a 
and b are the indices for different groups. There are many possible coexistent clustered states that depend on the 
initial conditions. With heterogeneity, clustered solutions of near synchronous neurons are possible as in the near 
synchronous solution. If the heterogeneity is nonuniform and clustered then groups may form more easily among 
neurons with similar intrinsic frequencies. 



C. Asynchrony, Harmonic Locking and Suppression 

Heterogeneity breaks the symmetry of the network. As we have seen above, if the heterogeneity is mild then phase- 
locked solutions can exist. For large enough heterogeneity there are no phase- locked solutions. When phase- locking 
does not exist several possibilities can arise. We have discussed asynchrony, harmonic locking and suppression for two 
coupled neurons. These states have their counterparts for a network of N neurons as well. 



If the heterogeneity is large then our analysis in Sec. VIA shows that stable synchrony can be broken. If the 



applied currents are very strong or the synapses are very weak then the period can become too short to support 
stable synchrony. In this case, the neurons will become effectively decoupled and fire asynchronously as in the two 
neuron case. This was observed numerically For heterogeneous neurons, the average density of neurons at a 

given phase need not be constant as in the homogeneous case. If the heterogeneity is very broad then it is possible 
that asynchrony and synchrony can coexist [p"8|jl9| ] . 

If the applied currents are weak and the synapses are strong then stable synchrony can again be broken. For 
inhibition, the possible outcomes are harmonic locking or suppression. Here the slower neurons cannot keep up with 
the faster ones. The dynamics will be similar to two self-inhibited neurons. The influence of the other neurons in the 
network with weakly heterogeneous applied currents acts like self-inhibition. If a neuron is too slow compared to a 
group of other neurons, then it can drop out of the rhythm. As the synaptic strength is increased, more and more 
neurons will drop out until only one remains. The condition for suppression can be estimated from the combined 
synaptic input the neuron receives over a period. Similarly, complicated harmonic patterns can form where the slower 
neurons lock to sub harmonics of faster neurons. In this case a large set of coexistent patterns exist |^l|. 



VII. DISCUSSION AND CONCLUSIONS 

We have examined the existence and stability of phased-lock solutions in a network of neurons with heterogeneous 
intrinsic frequency. With mild heterogeneity, we have shown that stable phase-locking is possible but it is fairly fragile. 
When periodic phase-locking breaks, states of asynchrony, harmonic locking or suppression may arise. In accordance 
with previous work ||l0|-[l^ , we have found that phase- locking tendencies are strongly determined by the type (I or 11) 
of the neuron model. We also distinguish between saturating and nonsaturating synaptic types and find that there 
can be differences in phase-locking properties. Neurons with saturating synapses tolerate the effects of heterogeneity 
better. 

Our analysis may clarify some of the apparent differences in previously derived conditions for synchronization with 
inhibition. We find that stable synchrony requires a nonzero rise time of the synaptic current as emphasized in 
Refs. |l^,^. This was implicitly assumed in Ref. jl2j. We also find that it requires that the rise time not be too slow 
which was emphasized in Ref. and assumed in Ref. |^^. However, Ref. also requires that the synaptic decay 
time not be too fast which was also observed in numerical simulations in Ref. |39[| . This condition may be implicit in 
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the behavior of the kernels. For some neuron models, the e kernel may become negative if the synaptic decay time 
becomes too fast. For example, if the threshold of the integrate-and-fire model changes with the synapse then this 
condition may be present . 

Our results rely on the spike response method to characterize the dynamics of a threshold crossing variable in terms 
of a set of response kernels. The existence and stability of phase-locked solutions are then established by the time 
courses of these kernels. Critical to our analysis is that the kernels exist and can be computed. For the integrate-and- 
fire model the kernels can be obtained explicitly. For more complicated models, they most likely must be obtained 
numerically p7| . However, we can still infer general properties on how they should behave. 

The neuron type (I or II) governs the general characteristics of the synaptic kernels. For example, positive inputs 
always advance the firing of Type I neurons. This implies that the synaptic kernel always remains positive. As noted 
previously l]!!],^, we find that Type I and Type II synapses will greatly affect the phase-locking tendencies of the 
neurons. Our formalism shows why this is the case from the shapes of the synaptic kernels. It should be noted that the 
expansion into kernels is not always guaranteed. In particular, the separation of a membrane kernel from a synaptic 
kernel may not always be possible if the neuron is strongly nonlinear. In such a case, an integral approximation can 
still to be obtained and linearizations around phase-locked solutions may be possible to examine stability. 

We emphasize that we have computed the local stability of the states. This is in contrast to previous strategies 
where the global stability is computed ^,^,|3j. In these works, the class of initial conditions that will approach a 
synchronous solution is found. An absorption criterion is used that assumes that once the oscillators are within a 
given phase of each other then they remain locked and no longer influence each other. In our formulation, the synaptic 
coupling can still influence locked neurons and be enough to destabilize the dynamics if certain other criteria are not 
met. In the case of excitatory coupling, we find that although oscillators may approach each other globally they 
may not be locally stable. It had been previously assumed that the reasons for disparate results for the stability of 
synchronized pulse-coupled integrate-and-fire oscillators were a result of the singular nature of the delta-distribution 
synaptic coupling [^jl2|. However, it is the assumptions regarding what constitutes a locked state that separates the 
conclusions. Ideally, a combined analysis should be done where a global stability analysis is used to examine the class 
of initial conditions that approach a phase-locked solution and then a local stability analysis is performed to see if 
the locked state is sustainable. 

We only considered heterogeneity in the applied currents which is a very simple form of disorder because it is 
not necessary to average over different instances of the heterogeneity. For instance if we choose N random current 
values from a fixed probability distribution and if N is large enough, then a single instance is equivalent to any other 
instance. This is because the network is insensitive to rearrangements of the current values between the neurons. We 
can always relabel the neurons in the network. Only the distribution of the gaps between current values are important 
and this is fixed for large enough N . This property allowed us to gain insight into the N neuron network from the 
two neuron network. On the other hand if the heterogeneity were to be included in the coupling strengths then a 
single instance of the heterogeneity may not give an accurate picture of what a typical random network will do. A 
rearrangement of the coupling strengths could greatly affect the network dynamics. Here, an understanding of two 
heterogeneously coupled neurons may not give any insight into the behavior of a large network. One must average 
over the heterogeneity to get an idea of what a typical network will do. Then one must probe the fluctuations around 
this average. For two neurons it seems that heterogeneity for synaptic strengths would have similar results to what we 
derived for heterogeneity in the applied current. Phase- locked solutions would be shifted away from the homogeneous 
solutions. Our stability results would then apply analogously. However, the behavior of large networks may be quite 
different from the behavior of two neurons. 
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APPENDIX A: CONDUCTANCE-BASED NEURON MODEL 

For numerical simulations we used a single-compartment model neuron with inhibitory synapses obeying first-order 
kinetics pl[ |. Membrane potential in each point neuron obeyed the current balance equation 

=h- iNa - Ik -II- Is, (A1) 

where C — l/iF/cm^, li is the applied current, /jva = gNa'm%:,h{Vi — Vjva) and Ik — gKn*{Vi — Vk) are the Hodgkin- 
Huxley type spike generating currents, II — gL{Vi — VL) is the leak current and Is — gsSj{t){Vi — Vs) is the synaptic 
current. The fixed parameters used were: g^va — 30 mS/cm^, — 20 mS/cm^, gj^ = 0.1 mS/cm^, V^va 45 mV, 
Yk = —75 mV, Vl = —60 mV, Ys — —75 mV. The phenomena described here seem largely independent of specific 
neuronal parameters. 

The activation variable m was assumed fast and substituted with its asymptotic value maoiv) — (1 + exp[— 0.08(^'-|- 
26)])~^. The inactivation variable h obeys 

dh _ feoo(f) - h 
at Th(v) 

with /loo(w) = (1 + exp[0.13(w + 38)])-\ Th{v) = 0.6/(1 + exp[-0.12(w + 67)]). The variable n obeys 

dn noo{v) — n 
dt Tn{v) 

with n^{v) = (1 + cxp[-0.045(w + lO)])-\ t„(z;) = 0.5 + 2.0/(1 + cxp[0.045(w - 50)]). 
The gating variable Sj for the synapse is assumed to obey first order kinetics of the form 

^^FiV,){l-S,)-(3S„ (A4) 

where F{Yj) = 1/(1 + cxp[— V,]), and Yj is the voltage of the pre-synaptic neuron and transmission delays are 
neglected. 

The ODEs were integrated using a fourth-order Runge-Kutta method. The free parameters were scanned across 
the following ranges: for applied current li, 0-10 /iA/cm^; for gs, 0-2 mS/cm^; for the synaptic decay time constant 
T,, 5-50 ms. 



(A3) 



APPENDIX B: PHASE-LOCKING OF TWO COUPLED NEURONS WITH NO MEMORY 

As an example, we consider a model where the response kernels have no memory beyond the previous firing. The 
phase-locked solutions are then given by 

f 1 (T + 0ir) = 1 = /i + 7?(T) + Je(T - (/.T), (Bl) 
V2iT + ^2T) = 1=I2 + v{T) + Je{<i>T) + Je{T + (j)T), , (B2) 
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This implies that G(0) = J[e{(t)T) + e(T + (j>T) - e(T - 0T)]. The perturbed dynamics Eq. are given by 

5{n) — M ■ d{n — 1), where 




r,{lT) 



Ji{lT-<pT) 

ill vi 

Je(lT+(l)T) Je(4,T)ri(lT) f){lT) i(4,T)e(lT-4>T) 

V2 V1V2 V2 



(B3) 



vi=T^{T) + Jk{T-c^T), 

V2 -= MT) + Je{<t>T) + Je{T + (f)T), 



(B4) 
(B5) 



and the prime refers to the first derivative with respect to time. By definition, v{t) must approach the threshold from 
below which implies that Vi > 0. 



The eigenvalues can be solved for exactly. We can verify that the rows of M{1) in (B3) sum to unity and so one 
eigenvalue is 1. The other eigenvalue is then A = TrM — 1 = det Af. For stability, we require |A| < 1 which implies 
that < TrM < 2 or -1 < det M < 1, where 



Tr M = 



J 



Vl 



V2 



V1V2 



(B6) 



detM 



viT) j^eiT - <j>T)eiT + ^T) 



V1V2 



V1V2 



(B7) 



We now apply our results of Sec. IV and compare to the exact eigenvalues. If all the elements of M are positive 
then we immediately see that stability is achieved. Since the rows must sum to 1 and the elements are positive then 
each element itself is less than 1. This immediately implies that < detM < 1 as required for stability. However, 
this condition although sufficient for stability is certainly not necessary. An examination of Tr M and det M shows 
that the elements of M can become negative and still maintain < Tr A/ < 2 required for stability. 

We next examine the necessary condition G'{4>) > of Theorem |^. We first note that G'{(f)) — JT[e{(f)T) + e'i{T + 
(j)T) + e(T - (f)T)]. Consider the synchronous solution for synchrony (0 = 0) where G"(0) = jr[e(0) + 2e(r)] > 
is clearly necessary for stability. If in addition, e(0) = and rj{T) > then an inspection of condition (B6) shows 
that G"(0) > is necessary and sufficient for stability. For the anti-synchronous solution {(j) = 0.5), G"(0.5) = 



JT[2e(0.5T) + e(1.5T)] > is necessary for stability from condition (B6) or (B7). However, unlike the synchronous 
case, it cannot be made into a sufficient condition without making additional nontrivial constraints on e(0.5r) and 
e{1.5T). 
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FIG. 1. Examples of possible states for two heterogeneous neurons from numerical simulations of the model given in Appendix 
a) Asynchrony with 7i = 9.0, I2 = 9.9 fiA/cm^; Qs = 0.25 mS/cm'^; Ts — 10 ms. b) Near-synchrony with Ji = 1.6, I2 = 
1.78 /iA/cm^; Qs = 0.25 mS/cm^; Ts — 10 ms. c) Harmonic locking with Ji = 9.0, I2 ~ 9.9 ^A/cm^; Qs = 0.5 mS/cm^; Ts = 
10 ms. d) Suppression with /i = 1.6, I2 ~ 1.78 /lA/cm^; Qs — 0.5 mS/cm^; Ts = 10 ms. 
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FIG. 2. Example kernels. An kernel with a spike is shown in a). The kernel for the integrate-and-fire model would be 
identical without the spike. In b) we show the e kernel from Eq. ( 2.7) for the saturating case (solid line) and the nonsaturating 
case (dashed line). Two examples of kernels for Type II neurons are shown in c). Note that they have both positive and 
negative parts. 
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